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Abstract 

We study a scalar hyperbolic partial differential equation with non-linear terms similar to 
those of the equations of general relativity. The equation has a number of non-trivial analytical 
solutions whose existence rely on a delicate balance between linear and non-linear terms. We 
formulate two classes of second-order accurate central-difference schemes, CFLN and MOL, for 
numerical integration of this equation. Solutions produced by the schemes converge to exact 
solutions at any fixed time t when numerical resolution is increased. However, in certain cases 
integration becomes asymptotically unstable when t is increased and resolution is kept fixed. 
This behavior is caused by subtle changes in the balance between linear and non-linear terms 
when the equation is discretized in space. Changes in the balance occur without violating a 
second-order accuracy of the discretization. We thus demonstrate that second-order accuracy 
and convergence at finite t does not guarantee a correct asymptotic behavior and long-term 
numerical stability. Accuracy and long-term stability of integration are greatly improved by an 
exponential transformation of the unknown variable. 
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1 Introduction 

It is well known that numerical integration of Einstein's equations in a 3+1 form often leads 
to instabilities and terminates prematurely. Some of the instabilities have been related to the 
presence of constraints, and some to the existence of gauge degrees of freedom in Einstein's 
equations (PP,|2], and references therein). These instabilities are intrinsic to the equations 
of general relativity (GR) themselves. Numerical instabilities may also arise due to a bad 
choice of a finite-difference scheme. 

In this paper we introduce a scalar hyperbolic partial differential equation with non-linear 
terms similar to those of much more complicated equations of GR (Section 2). The equation 
has a number of non-trivial analytical solutions (Sections 3 and 4). We experiment with 
several finite-difference schemes for numerical integration of this equation (Section 5). 

Solutions produced by the schemes converge to exact solutions at any fixed moment of 
time t when numerical resolution is increased. However, in certain cases numerical integra- 
tion becomes unstable when t is increased while the resolution is kept fixed (Section 6). We 
trace this behavior to a special structure of the non-linear terms and demonstrate that dis- 
cretization of non-linear terms leads to a finite-difference system whose asymptotic behavior 
may differ qualitatively from the corresponding behavior of a continuous system (Section 
7). The accuracy and stability of integration can be greatly improved by an exponential 
transformation of the unknown variable (Section 7). 

2 A model equation 

Consider a scalar, quasi-linear, hyperbolic partial differential equation of two independent 
variables, t and x, 



where g = g(t, x) is the unknown, and constant coefficients, 011022 — 021 < 0- The 

non-linear term in (|2.1|) has the form g^ 1 ^ (first derivatives squared). 

The reason for considering ()2.1|) becomes obvious if we recall that the structure of a Ricci 
tensor R a ^ is R ~ ^ dT + ^ IT, where Cristoffel symbols T ~ g~ 1 dg, and g is the metric. 
Thus, Rob is represented as a sum of the terms R ~ g~ 1 d 2 g + g~ 2 (dg) 2 . Equation (|2.1j) 
thus mimics a type of non-linearity present in GR equations R a b = 0. In particular, (|2.1j) 
may have a singularity related to g becoming zero. 

Equation (|2.1|) can be reduced to its normal form, 



by a linear transformation which preserves a quadratic form of the non-linearity. We work 
below with a simpler equation (J2.2j) . This is a non-linear hyperbolic equation with a char- 
acteristic speed equal to 1. 

By introducing a new variable K = g t , we rewrite ()2.2|) as a system of two first-order in 
time partial differential equations 




(2.1) 





9t = K, 
Kt = 9xx ~ 



g-\aK 2 + f3g 2 x + 1 Kg x ) 



(2.3) 
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which resembles the evolutionary part of GR equations in a standard ADM 3+1 form (with 
zero shift and constant lapse). By introducing yet another variable D = g x we further rewrite 
system of first-order PDE 

9t = K, 

K t -D x = R, (2.4) 
A - K x = 0, 

where 

R= -g-\aK 2 + (3D 2 + jDK). (2.5) 
This system has a complete set of real eigenvalues and eigenvectors. 

3 Exponential transformation 

A transformation 

9 = et (3.1) 

where is a new unknown, will play an important role in subsequent sections. After substi- 
tution of (j3.1|) into (|2.2j) we obtain a PDE for a new unknown 0, 

<fa = <f>** - (« + 1)0? - (0 - 1)0' - 70x0*- (3.2) 

The transformation (J3.1)) removes g~ x multiplier in front of the non-linear term in (J2.2j) . and 
maps < g < oo onto — oo < (j) < oo so that values of g < are excluded. This is consistent 
with GR where three-dimensional metric of space-like hypersurfaces is positive-definite. 
We can introduce new variables 

iff = <jh, = <fr x , (3.3) 
and rewrite ()3.2|) in the equivalent two-equation form similar to (|2.3J) as 

i2 , , (3-4) 



<pt = 1p, 




i>t = <frxx - (a + l)tp 2 - 


-(P- 


form similar to ()2.4|) as 




<f>t 


= v, 


A - o x 


= s, 


t - 


= 0, 



(3.5) 



where 

S = _( a + i)^, 2 - (J3 -l)$ 2 -^9. (3.6) 
System ()3.5|) obviously has the same eigenvectors and eigenvalues as (|2.4j) . 
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4 Analytic solutions 

For a choice of parameters 

a = -l, (3 = 1, and 7 = 0, (4.1) 
()3.2|) becomes a linear hyperbolic PDE 



>tt — Yxx 



(4.2) 



whose general solution is 

<P = <f) 1 {x + t)+<f) 2 {x-t), (4.3) 
where 0^2 are arbitrary functions. The original equation (J2.2j) remains non-linear, 

gtt-g xx = g~\g 2 t -gl). (4.4) 

Its general solution is 

g = gi {x + t) -g 2 (x-t), (4.5) 

where g±, g 2 are arbitrary functions. 

For a, /3, and 7 other than (|4.1j) . equation (J3.2j) is non-linear but we can find particular 
solutions of this equation in a form of a wave running with a constant speed a, 

= 0(0, ( = x + at. (4.6) 
Substituting ()4.6|) into ()3.2|) we obtain an ordinary differential equation 

where 

6= (a + l)a 2 + /3-l + 7a. (4.8) 

We need to solve this equation for <fi. 

We consider separately the cases b(a 2 — 1) 7^ 0, b = 0, and a 2 = 1. If 6(a 2 — 1) ^ 0, 
integration of ()4.7|) gives 

a 2 - 1 

4> = p ln(c + x + at) + ci, where p= — - — , (4.9) 



and c and C\ are arbitrary constants. The corresponding solution of ()2.2|) is 

g = e Cl (c + x + at) p . (4.10) 
If a 2 7^ 1 and 6 = 0, which will happen if we chose 



- 7± ^+4(q+ !)(!-£) 

2(a+l) r V ' 
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then integration of ()4.7|) gives 

<f> = Co(x + at) + ci, (4.12) 
and the corresponding solutions of ()2.2j) are 

g = e co(x+at)+ C1 _ ^jgj 

Finally, if a 2 = 1, integration of ()4.7|) gives 

0i(x + at) if 6 = a; + /5 + 07 = 1 a 

const if6 = a + /? + a7^0 ' [ ' 

where <p\ is arbitrary. The speed a of solutions (j4.10j) . ()4.13|) differs from the characteristic 
speed of (j2.2j) . These running wave solutions exist as a result of a delicate balance of linear 
and non-linear terms in ()2.2|) . 



5 Numerical schemes 

Let us now consider how ()2.2|) may be solved numerically. Without non-linear terms, ()2.2|) 
is a scalar wave equation g tt = g xx . A classical, central-difference, second-order accurate 
scheme for this equation was introduced in 1928 by Courant, Friedrichs and Levy (see |3] 
and chapter 10 in ji]), 

At 2 Ax 2 ' 1 ' 

where g™ are determined at mesh points Xi = iAx and t n = nAt. The scheme is stable for 
Courant numbers cfl — ^ < 1. 

We now proceed to expand (|5.1|) to the non-linear equation. We cast (jO| into a numer- 
ically equivalent two-equation form, 

KT l -KT h _ g? +l -2g? + gl 1 

At ~ Ax 2 ' (5.2) 

„n+l „n -, 
9j ~ 9j K n +\ 

At 1 ' 

and extend it to (J2.3)) by adding a discretized non-linear term to the right-hand side of the first 
equation in (|5.2j) . In order to maintain an overall second-order accuracy of the scheme, we 
must add the non-linear term evaluated with second-order accuracy at grid points (xi,t n ). 
This requires second-order accurate values of g and K at these points. Values of g™ are 

already defined there but are not. Taking = |(i^™ +2 + 2 ) will give us the 
desired accuracy but will also render the scheme implicit. To escape this difficulty, we use a 
predictor-corrector approach. 

We write a discretized non-linear term (|2.5jl as 

K{g h K u A) = -9i l {aK? + (3D 2 + 1 K l D l ) , (5.3) 
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where Di = z£x(9i+i — 9i-i), without explicitly specifying a moment of time at which ^ and 
Ki, must be taken. Using this notation, we write a second-order accurate explicit predictor- 
corrector scheme for ()2.3|) as 



CFLN1 



Ki = K 



At 



i+l 



2g? + g 



i-l 



Ax 2 



K7 : 
n+1 
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= K t + ^ (n(g?, K h A n ) - ng?, K~\d? 

g? + AtK" +1 >, 



n(g^K t \Df) J (predictor) 
(corrector), 



(5.4) 

We will refer to (|5.4j) as to a CFLN1 scheme. Initial conditions for the scheme must be 
provided at t = t° for g>j and at t = t° — |At for K { . Boundary conditions are required 
only for g and must be provided at t n . We show elsewhere that the scheme is stable for 
cfl < 1 according to a standard VonNeuman stability analysis J5j. Thus, it may be expected 
to converge at any t to exact solutions when the resolution is increased, Ax — > 0. Numerical 
experiments presented in the next section confirm this assertion. 

Note, that CFL1 can be easily cast into a three-equation form, a numerical counterpart 
of (|231) 



CFLN2 : 



K; - K 



D r 



i+7 



D n 



At 



Ax 



TL^Xi \m (predictor), 



K n +2 



Ki + -At [ K(g?, K h D?) - K(g?, K, 



(corrector), 



D 



n+l 



K 



i+l 



K 



At 



Ax 



-,n+l 



At 



k: + \ 



where we introduced new quantities 



9i+i ~ 9i 
Ax 



(5.5) 



(5.6) 



so that Di = ^(D i+ i + D^i). We will refer to ()5.5j) as to a CFLN2 scheme. Schemes 
CFLN1 and CFLN2 are numerically equivalent and have identical stability and convergence 
properties. 

A second group of numerical schemes considered in this paper is based on a method-of- 
lines approach. We discretize f)2.3j) in space using central differences, 



MOLlfn) 



d Jl= K 
dt 17 
OK 



dt 



9i+i +9i-i- 29i 
Ax 2 



+ K{g u K^ Di 



(5.7) 



and integrate the resulting system of ordinary differential equations in time using Runge- 
Kutta methods of orders n = 2,3,4. In what follows, we refer to (j5.7J) as to an MOLl(n) 
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scheme. The scheme is second-order accurate and stable for sufficiently small Courant num- 
bers. Initial conditions for g^ and Ki must be provided at t — t°. Boundary conditions are 
required for gi at times t n . 

Instead of (|2.3J1 . we can discretize (|2.4j) in space, 



M0L2 : 



r 9_g 1 = 

at 11 



i+1 



dt Ax 
OK D i+ i - D. i 



(5.8) 



dt 



Ax 



+ n(gi,Ki,Di), 



and then use Runge-Kutta methods of order n to integrate ()5.8|) in time. We will refer 
to O as to a M0L2(n) scheme. Using (CTRjl . it is easy to verify that M0L1 and M0L2 
schemes are in fact equivalent. 

Instead of a Runge-Kutta, one can use any other stable ODE integrator in M0L1 and 
M0L2, e.g., an iterative Crank-Nicholson (ICN) scheme with appropriate number of itera- 
tions 6J. We do not present here results obtained with the ICN as they are similar to those 
obtained with both CFLN1 and Runge-Kutta M0L1 schemes, and do not change any of the 
conclusions of the paper. 



6 Numerical convergence and asymptotic stability 

In this section, we present examples of numerical integration of ()2.3|) . They illustrate conver- 
gence of the schemes at fixed time t when the resolution is increased, Ax — > 0. The examples 
also illustrate numerical difficulties which may arise when Ax is kept fixed and integration 
time is increased, t — > oo. 

As a first example, consider equation (|2.2|) with a = — |, /3 = |, 7 = 0. We pick a wave 
speed a = 2 and find from (|4.8|) . (|4.9|) . (|4.1(J|) a particular growing solution 

El: g = (x + 2t)^. (6.1) 
For a = 0.1 we find a decaying solution 

/ t \ " 3 - 88 

E2 : g=[x + -j . (6.2) 

We wish to obtain solutions El and E2 numerically on interval 0.1 < x < 1.1, and t > 
using iV grid points with coordinates 

Xi = 0.1 + Ax (i- 1/2) , Ax = l/N, (6.3) 

and boundary conditions 

A t A t 

g% =0(0.1-— ,nAt), g n N+ ,=g(l.l + —,nAt), (6.4) 
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where g is the corresponding exact solution (El or E2). Results of numerical integration are 
presented in Tables 1 and 2, and in Figure 1. 

Table 1 shows convergence of a MOLl(4) scheme for a growing solution El. The maximum 
error norm L ro indicates a second-order convergence. Relative error of integration decreases 
with time. It is possible to continue stable integration of El until the limit of large numbers 
is reached in a computer. Results for M0L1(2), M0L1(3), and CFLN1 schemes are similar. 

Table 2 illustrates convergence of a M0L1(4) scheme for a decaying solution E2. Similar 
to El, convergence is second-order. However, the relative error is much larger than in the 
El case and it grows with time. 



N 


^oo(^l) 




16 


2.4E-04 


7.9E-05 


32 


4.7E-05 


1.4E-05 


64 


1.2E-05 


5.9E-06 


128 


3.0E-06 


1.6E-06 


256 


7.5E-07 


4.2E-07 



Table 1: Convergence of MOLl(4) numerical scheme for a growing solution El. Integration is carried 
out with a time step At = \Ax for 0.1 < x < 1.1 and t > 0. The error norm = maxj \gi/g e — 1 1 , 
where g e is the exact solution ()6.1j) . is given for two moments of time, t\ = 9.9 and t 2 = 24.75. 

Figure 1 compares numerical solutions in the middle of the interval, x = 0.6, to the exact 
solution E2 at this point. For t < 20, the N = 2048 numerical solution and the exact solution 
cannot be distinguished on the plot. The iV = 128 numerical solution shows large deviations 
thich grow with time. For N = 64, the deviations are so violent that the code crashes well 
before reaching t = 20. Results obtained using MOLl(2), MOLl(3), and CFLN1 are similar. 
By increasing N we can achieve a convergent solution at any fixed time t. However, if we 
keep the resolution fixed and increase t, we find that numerical errors make a long-term 
stable integration impossible. 



N 


-^oo(^l) 


Loo{t2) 


64 


2.9E-01 


NaN 


128 


1.3E-01 


2.2E-01 


256 


3.0E-02 


8.1E-02 


512 


7.5E-03 


1.9E-02 


1024 


1.9E-03 


4.8E-03 


2048 


4.7E-04 


1.2E-03 



Table 2: Convergence of MOLl(4) numerical scheme for a decaying solution E2. Integration is 
carried out with a time step At = \Ax for 0.1 < x < 1.1 and t > 0. The error norm L ro = 
maxj \gi/g e — 1|, where g e is the exact solution Q6.1|) . is given for two moments of time, t\ = 9.9 
and t 2 = 24.75. 
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t 



Figure 1: A decaying solution Q6.2|) at x = 0.6 as a function of time. Obtained using MOLl(4) 
scheme and cfl = 1/2. Solid lines - numerical solutions for N = 128 and N = 2048. The exact and 
N = 2048 numerical solutions cannot be distinguished on this plot. 
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N 


Loo(tV) 


L 00 (t2) 


16 


2.6421E-04 


2.2272E-03 


32 


5.6685E-05 


4.4543E-04 


64 


1.2938E-05 


2.0090E-04 


128 


3.0747E-06 


5.9922E-05 


256 


7.3409E-07 


1.5723E-05 



Table 3: Convergence of numerical solutions g_ for two moments of time t\ = 1.2375 and i 2 = 
3.7125. 



To understand this behavior of numerical schemes, it is instructive to consider a special 
case of exponential solutions ()4.13|) which can be investigated analytically. As an example, 
we take a = 1, /3 = 7 = in ()2.3|) . This choice of parameters gives a pair of exponential 
solutions ()4.13|) with the speed a = ±775 > 

g± = exp ( x ± —=) . (6.5) 



We now attempt to obtain g± numerically on the interval < x < 1 discretized using N grid 
points, 

x i=T?( i -lh i = l,...,N, (6.6) 



N \ 2, 

and applying Robin boundary conditions = 1 for i = and i — N + 1, 

5-0 = exp(hi5( 2 - 2Axlii0i), t/jv+i = exp(ln#jv-i + 2Ax\ng N ). (6.7) 

We use these boundary conditions because, in this particular case, they completely eliminate 
the influence of boundaries on a numerical solution in internal points (see below). For 
integration in time we use MOLl(4) with cfl = |. By varying Courant number in the 
range 1 < cfl < and Runge-Kutta order from n = 2 to n = 4 we verified that errors of 
integration in time in our numerical experiments are less than 10~ 4 of the truncation errors 
introduced by the spatial discretization (|5.7j) . 

Results of numerical integration for a decaying solution g_ are shown in Figure 2 for 
resolutions iV = 16 through iV = 1024. In the beginning, numerical solutions follow the 
exact solution but eventually begin to deviate and grow exponentially. With increasing N, 
the exponential growth starts later. However, time t s of stable integration is proportional 
only to a logarithm of a number of grid points, t s ~ IniV. Integration beyond, say, t ~ 10 
would require an unpractical larger number iV > 10 6 . Table 3 illustrates convergence of 
numerical solutions for two different moments of time, t\ = 1.2375 and £2 = 3.7125, which 
are both less than t s . There is a second-order convergence at these times ( integration of g + 
does not present any difficulties and can be continued indefinitely). 

We now analyze the reason for an asymptotic instability discussed above. Using initial 
conditions #i(0) = exp(xi), we can present numerical solutions gt(t) at t > as 

dfi 

9i(t) = fi{t)exp(xi), Ki{t) = -^exp(xi), (6.8) 
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Figure 2: A decaying solution g_ as a function of time. Obtained using a MOLl(4) integrator and 
cfl = 1/2. Solid lines - numerical solutions with N = 16 through N = 1024. Dashed line - exact 
solution. 
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where fa are functions of time. Substituting ()6.8|) into ()5.7|) . we obtain a set of identical 
equations for all fi, 

^ = c/,-if- (6.9) 

where 

Ax I „-Ai _ o 

c =- — If- - = l + 0(Ax 2 )>l (6.10) 

Aar 

is a constant independent of i (Robin boundary conditions ()6.7|) are necessary to ensure that 

(16. 9 j) holds for the outermost grid points i = 1 and i = N). Initial values fi{0) = 1 and 
dfdo) 
at 

ordinary differential equation 



;) , - are also identical for all i. We thus can ignore index i in (|6.9)1 . and use a single 

to describe numerical solutions = /(£) • ^(O) at all grid points. Integration of (|6.11|) 
gives 



df\ _cf 2 C 



+ ( 6 - 12 ) 




where C = const. We set it to C = < to satisfy initial conditions, and finally obtain 
the following differential equation 

At \ „ti 1 

(6.13) 

The value of c = 1 corresponds to a continuum limit of Ax — > 0. Difference in the behavior of 
numerical and analytical solutions comes from the presence of the term . Note, that this 
term is 0(Ax 2 ). Its variations do not change the second-order accuracy of the algorithm. 
The term is negative since c > 1. 

Consider first an exponentially growing solution g + . For this solution, > 0. There- 
fore, we must initially take a + sign in (|6.13|) . The expression under the square root in (|6.13J) 
is positive at t = and it will only increase with increasing / because the term — > 
when / — > oo. The relative difference between numerical and analytical solutions will tend 
to zero when t — > oo as well. 

For a decaying solution <?_, asymptotic behavior of exact and numerical solutions is 
qualitatively different. For the exact solution, we have K — > and / —>■ when t —>■ oo. 
But from ()6.12|) we observe that for c other than zero / cannot become zero because 
in ()6.12|) has a minimum |£ = at a finite / = f s — (2(c — l)/^ 1 / 4 > 0. For a decaying 

solution, a y < and we must initially take — sign in (|6.13p . When / reaches the value of 
f s we have df /dt — but the second derivative remains positive 

ff- =2c f + (c-l)/f>0. (6.14) 
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As a result, the numerical solution switches at / = f s from — branch to + branch in ([6. 13)1 . 
and begins to increase, / — ► oo when t — > oo. 

We can estimate the moment of time, t s , when the solution switches from the — to the 
+ branch from 

e~% ~ f s = (2( c - l)/c) 1/4 ~ Ax 1/2 . (6.15) 

When the resolution is increased, f s becomes smaller and is reached at later time. But 
(jSHH) shows that t s oc (In Ax) -1 . To increase a period of stable integration t s for a decaying 
solution, we must decrease Ax (increase N) exponentially! 

If we impose boundary conditions other than 1)6. 7J) . deviations from the exact solution in 
internal points grow as described by (J6.13)) only until a signal from the boundary reaches 
them. After that, interaction with the boundary leads to a violent instability and a termina- 
tion of numerical calculations. The reason for an instability observed for a decaying solution 
E2 (Figure 1) appears to be the same. Truncation errors lead to an imperfect balance of 
linear and non-linear terms, and to a deviation of a numerical solution from the exact one. 
Subsequent interaction with the boundaries amplifies these deviations, and the calculation 
eventually terminates. 



7 Integration in logarithmic variables 



We now show that logarithmic variables <p — lag, if> — <ft t , and 9 = ip x allow to significantly 
improve the accuracy and stability of numerical integration. 

Equation ()3.2|) for can be integrated numerically using the same schemes CFLN1 and 
MOLl(n). We must simply replace g, K and D in these schemes with 0, if), and 9, and 
to substitute the non-linear term TZ with its logarithmic counterpart S (|3.6|) . The schemes 
than become 



CFLN1 : 



and 



where 



and 



- 20F + 



'i-i 



Ax 2 



(predictor), 



At 



MOLl(n) 



^ + T (s(^,a-w; 

C + At< + ^, 



dt 

9 f .- 



9™)\ (corrector), 



Ax 2 



+ S(i) i ,9 i 



2Ax 



S i = -(a + l)tf-(j3-l)%- T (>i0i- 
Schemes CFLN2 and MOL2 can be rewritten in a similar way. 



(7.1) 
(7.2) 

(7.3) 
(7.4) 
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5 



10 



15 



20 



25 



Figure 3: Decaying solutions 1)6. 2J) at x = 0.6 as a function of time obtained in logarithmic variables 
using CFLN1 scheme with cfl = 1. Solid lines - numerical solutions for N = 128 and N = 2048. 
Dashed line - exact solution. The solutions cannot be distinguished on the plot. 

Let us first consider how a CFLN1 scheme will reproduce an exponentially decaying 

solution g_ (J6.5j) . For this solution, initial conditions linear function of x, and 

_ i 

ipi 2 = a constant. With these conditions, the right-hand sides of the first two equations 
in (j7.1j) become zero, and it is easy to verify that ipi remain constant through all subsequent 
time steps, whereas <pi decrease with time linearly, 0™ = Xi — z ^=. The same is obviously true 
for MOLL Transformed to logarithmic variables, CFLN1 and MOL1 reproduce exponential 
solutions exactly. 

Next, we test the reformulated schemes on a decaying solution (J6.2j) . Figure 3 shows 
numerical solutions at x — 0.6 as a function of t obtained using CFL1 scheme, and it 
must be compared to Figure 1. Table 4 illustrates convergence of reformulated CFLN1 and 
MOL1 schemes with increasing iV and should be compared with Table 2. We see from the 
comparison that logarithmic variables greatly reduce errors of numerical integration. 
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N 


Loo(tl) CFLN1 


L 00 (*2) CFLN1 


Loo(tl) MOL1 


Loo^) MOL1 


32 


3.9E-02 


8.9E-02 


4.1E-02 


2.1E-01 


64 


1.3E-02 


3.2E-02 


1.3E-02 


3.4E-02 


128 


3.4E-03 


1.1E-02 


3.7E-03 


1.3E-02 


256 


8.7E-04 


3.4E-03 


9.1E-04 


3.6E-03 


512 


2.1E-04 


8.9E-04 


2.1E-04 


9.0E-04 


1024 


5.2E-05 


2.2E-04 


5.2E-05 


2.2E-04 



Table 4: Convergence of numerical solutions 1)6.2(1 for two moments of time t\ = 20 and ti = 50 
using logarithmic variables and two numerical schemes, CFLN1 with CFL number cfl = 1 and 
MOLl(4) with cfl = 0.5. The error norm Loo — rn ax i|<?i/ 1 9e — 1|> where g e is the exact solution. 



8 Conclusions 

In this paper we introduced a scalar wave equation with non-linear terms similar to those 
of more complex equations of general relativity. The equation has a number of non-trivial 
analytical solutions useful for testing numerical schemes. 

We formulated two classes of finite-difference schemes for numerical integration of this 
equation. One (CFLN) is a non-linear extension of a classical second-order central difference 
scheme for a linear scalar wave equation j^j. Another (MOL) is based on a method-of-lines 
approach. The schemes have a comparable accuracy but MOL requires a larger number of 
right-hand side evaluations. 

Both schemes converge to exact solutions at any fixed t when numerical resolution is 
increased, Ax — ► 0. For some of the solutions, however, integration becomes unstable when 
resolution is kept fixed and t is increased. We trace this behavior to deviations from a 
perfect balance between linear and non-linear terms caused by discretization. As a result, 
the asymptotic behavior of numerical solutions may differ qualitatively from the asymptotic 
behavior of the corresponding exact solutions of a partial differential equation. An important 
point is that these deviations happen without violating a second-order accuracy. Therefore, 
having both convergence at finite t and an asymptotic instability is not a contradiction. 

An asymptotic instability seems not to be a fault of a particular numerical scheme. All 
schemes tested in this paper display this phenomenon regardless of how numerical solutions 
are advanced in time. CFLN, Runge-Kutta, and ICN-type integration results in the same 
asymptotic instability. We have no reason to believe that this phenomenon should be limited 
only to scalar non-linear equations. Examples presented in the paper clearly demonstrate 
that second-order accuracy of spatial discretization, although necessary for obtaining con- 
vergent numerical solutions at finite t, is not a guarantee of a correct asymptotic behavior 
of a numerical scheme. 

Finally, we have shown that an exponential transformation (|3.1|) leads to a significant 
improvement in accuracy and stability of numerical algorithms. We do not know if some 
of the difficulties encountered in numerical general relativity steam from a similar non- 
linear asymptotic instability, and whether integration of GR equations can be improved 
by an exponential transformation of variables. We believe that these questions are worth 
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investigating. In the accompanying paper we show how an exponential transformation of 
variables can be carried out for tensorial equations of GR. Numerical schemes formulated 
in this paper can be extended to GR equations written in both second-order (CFLN1 and 
MOL1) and first-order forms (CFLN2 and MOL2). 
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